#ifndef _CCD_H_
#define _CCD_H_
#include"system_hf.h"
#include<gsl/gsl_sf_coupling.h>
#include<vector>

#include<eigen3/Eigen/Eigen>
using Eigen::MatrixXd;



struct ccd_channel
{
  ccd_channel(){}
  ccd_channel(int _Tz,int _Par,int _J):Tz(_Tz),Par(_Par),J(_J){}
  int Tz,Par,J;
  vector<int> bra;
  vector<int> ket;
  MatrixXd Mat;
};

struct position//the channel and indice of state |a,b,J>
{
  vector<int> channel;//use vector to store different J with same a,b, the same below 
  vector<int> indice;
};

class CCD
{
 public:
  typedef System_Tz_HF System_HF;
  typedef System_HF::OrbitalType HFOrbital;
  typedef vector<HFOrbital> HFOrbitals;
  typedef TwoBodyState_TzParJ TwoBodyStateType;
  
 CCD(const System_HF* Ptr):pSystem(Ptr),HFOrbs(Ptr->Orbitals){setup();}

  double cross_coupled1(int a,int b,int c,int d,int J) const;
  
  void generateTwoBodyStates();
  void generateBlockMatrices();
  void setup();
  bool Ladders();
  bool Full();

  double gethhpp(int i,int j,int a,int b,int J,const vector<ccd_channel> & hhpp);
  double getpppp(int a,int b,int c,int d,int J,const vector<ccd_channel> & pppp);
  double gethhhh(int i,int j,int k,int l,int J,const vector<ccd_channel> & hhhh);
  vector<TwoBodyStateType> TwoBodyStates_hh;
  vector<TwoBodyStateType> TwoBodyStates_pp;
  vector<TwoBodyStateType> TwoBodyStates_phm;
  vector<TwoBodyStateType> TwoBodyStates_hpm;
  Matrix<position,Dynamic,Dynamic> CCD_position;
  vector<ccd_channel> CCD_V_hhhh, CCD_V_hhpp,CCD_V_pppp,CCD_T_hhpp,CCD_e_hhpp,CCD_V_hpmhpm,CCD_V_phmhpm;

  double TwoOrder_2p2h();

  
  double CCD_deltaE;
  const System_HF * pSystem;
  const HFOrbitals & HFOrbs;
  int FermiSurface;
  int totalOrbitals; 
};
#endif
